#include <stdexcept>
#include "../Okada.hxx"
#include "Iteration_Params.hxx"

Okada::Displacement Okada::source_contribution
(const bool &real_source,
 const double xi[],
 const FTensor::Tensor1<double,3> &coord)
{
  /* Real-source contribution */

  FTensor::Index<'a',3> a;
  FTensor::Index<'b',3> b;
  FTensor::Index<'c',3> c;

  FTensor::Tensor1<double,2>
    offset_coord(coord(1),origin(2) + (real_source ? coord(2) : -coord(2)));

  FTensor::Tensor1<double,2> transformed_coord;
  transformed_coord(0)=offset_coord(0)*z_transform_dip(1,1)
    - offset_coord(1)*z_transform_dip(1,2);
  transformed_coord(1)=-offset_coord(0)*z_transform_dip(1,2)
    + offset_coord(1)*z_transform_dip(2,2);

  const double p(transformed_coord(0));
  const double q(clamp(transformed_coord(1)));

  double eta[2]={clamp(p+width),clamp(p)};

  /* Reject singular cases */
  /* On fault edge */

  if(q==0 && ((xi[0]*xi[1]<=0 && eta[0]*eta[1]==0)
              || (eta[0]*eta[1]<=0 && xi[0]*xi[1]==0)))
    throw std::runtime_error("Bad input: On fault edge");

  /* On negative extension of fault edge */

  double r12(sqrt(xi[0]*xi[0] + eta[1]*eta[1] + q*q));
  double r21(sqrt(xi[1]*xi[1] + eta[0]*eta[0] + q*q));
  double r22(sqrt(xi[1]*xi[1] + eta[1]*eta[1] + q*q));

  bool kxi[]={xi[0]<0 && r21+xi[1]<eps, xi[0]<0 && r22+xi[1]<eps};
  bool keta[]={eta[0]<0 && r12+eta[1]<eps, eta[0]<0 && r22+eta[1]<eps};

  Displacement result;
  result.first(a)=0;
  result.second(a,b)=0;

  for(size_t k=0;k<2;++k)
    for(size_t j=0;j<2;++j)
      {
        const Iteration_Params p(xi[j],eta[k],q,kxi[k],keta[j],
                                 rotation_dip(2,1),rotation_dip(2,2));
        Okada::Displacement u_a(p.ua(dislocation,alp1,alp2,rotation_dip(2,1),
                                     rotation_dip(1,1))), u_temp;

        FTensor::Number<2> two;
        if(real_source)
          {
            u_temp.first(a)=-rotation_dip(a,b)*u_a.first(b);
            u_temp.second(a,b)=-u_a.second(a,c)*rotation_dip(b,c);
            u_temp.second(two,a)*=-1;
          }
        else
          {
            Okada::Displacement u_b(p.ub(dislocation,alp3,
                                         rotation_dip(2,1),
                                         rotation_dip(1,1),
                                         cos_dip2,sin_dip2,sin_cos_dip)),
              u_c(p.uc(coord(2),dislocation,alp4,alp5,rotation_dip(2,1),
                       rotation_dip(1,1),cos_dip2,sin_dip2,sin_cos_dip));
            u_temp.first(a)=rotation_dip(a,c)*(u_a.first(c)+u_b.first(c))
              +coord(2)*z_transform_dip(a,c)*u_c.first(c);

            u_temp.second(a,b)=
              rotation_dip(b,c)*(u_a.second(a,c)+u_b.second(a,c))
              +coord(2)*z_transform_dip(b,c)*u_c.second(a,c);
            u_temp.second(two,a)+=z_transform_dip(a,b)*u_c.first(b);
          }
        int sign=(j==k ? 1 : -1);
        result.first(a)+=sign*u_temp.first(a);
        result.second(a,b)+=sign*u_temp.second(a,b);
      }
  return result;
}

